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Abstract. The dynamical interactions of planetary systems may be a clue to their formation 
histories. Therefore, the distribution of these interactions provides important constraints on 
models of planet formation. We focus on each system's apsidal motion and proximity to dynam- 
ical instability. Although only ^^25 multiple planet systems have been discovered to date, our 
analyses in these terms have revealed several important features of planetary interactions. 1) 
Many systems interact such that they are near the boundary between stability and instability. 
2) Planets tend to form such that at least one planet's eccentricity periodically drops to near 
zero. 3) Mean-motion resonant pairs would be unstable if not for the resonance. 4) Scattering 
of approximately equal mass planets is unlikely to produce the observed distribution of apsidal 
behavior. 5) Resonant interactions may be identified through calculating a system's proximity to 
instability, regardless of knowledge of angles such as mean longitude and longitude of periastron 
(e.g. GJ 317 b and c are probably in a 4:1 resonance). These properties of planetary systems 
have been identified through calculation of two parameters that describe the interaction. The 
apsidal interaction can be quantified by determining how close a planet is to an apsidal sepa- 
ratrix (a boundary between qualitatively different types of apsidal oscillations, e.g. libration or 
circulation of the major axes). This value can be calculated through short numerical integra- 
tions. The proximity to instability can be measured by comparing the observed orbital elements 
to an analytic boundary that describes a type of stability known as Hill stability. We have set up 
a website dedicated to presenting the most up-to-date information on dynamical interactions: 
jhttp:/ / www.lpl.arizona.edu/~rory / research /xsp /dynamics] 



1. Introduction 

One of the most striking differences between known exoplanets and the giant planets 
of our Solar System involves the observed orbits: Observed exoplanets tend to have large 
eccentricities e, and small semi-major axes a, whereas the gas giants of the Solar Sys- 
tem have small e and large a. Recently, however, it has been shown that the dynamical 
interactions in many multiple planet systems (including the Solar System) show certain 
features in common (Barnes & Quinn 2004; Barnes & Greenberg 2006a [BG06a]; Barnes 
& Greenberg 2006c [BG06c]; Barnes & Greenberg 2007b [BG07b]). These shared traits 
suggest that the character of dynamical interactions (over 10"^ - 10^ years), rather than 
the present orbits, may be a more meaningful constraint on the origins of planetary sys- 
tems (Barnes & Greenberg 2007a [BG07a], BG07b). Considerations of shared dynamical 
properties have even resulted in the first successful prediction of the mass and orbit of 
an extrasolar planet, HD 74156 d (predicted by Barnes & Raymond [2004] and Raymond 
& Barnes [2005]; found by Bean et a/.2008). 

Now with over 200 extra-solar planets known, including ^ 25 multi-planet systems, 
we can look at this population as a whole, with increasing confidence regarding which 
common characteristics may be more than statistical flukes. We have identified com- 
monalities among planetary interactions that may be key constraints on the origins of 
planetary systems. We have identified parameters that quantify an interaction's prox- 
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imity to boundaries between qualitatively different types of motion. These two types 
of boundaries are the "apsidal separatrix" and the dynamical stability boundary. The 
apsidal separatrix is the boundary between different types of apsidal oscillations, e.g. 
libration and circulation. The stability boundary separates regions in which all planets 
are bound to the host star from those in which at least one planet is liable to be ejected. 
These two parameters are fixed quantities that do not vary over time, but they constrain 
how the systems evolve. Systems tend to lie close to these two boundaries. 

In this chapter we review our derivation of e, which quantifies proximity to an apsi- 
dal separatrix, in § 2. Then we describe (3, which parameterizes a system's proximity to 
dynamical instability, in § 3. In § 4 we present the observed distributions of these quan- 
tities and use them to constrain formation models. Finally in § 5 we draw our general 
conclusions. 

2. The Apsidal Separatrix 

Eccentricities in multiple planet systems oscillate due to secular interactions (see e.g. 
Laughlin & Adams 1999; Rivera & Lissauer 2000; Stepinski et aZ.2000; Michtchenko & 
Malhotra 2002; BG06a). Thus a currently observed eccentricity may not be its average 
value over the ^^10,000 year secular period. Validation of formation models through the 
eccentricity distribution is therefore inadequate; the eccentricity oscillations are a better 
description of a multiple-planet system's properties. 

We consider analytic and numerical models of planet-planet interactions to describe 
the oscillations. The analytic approach is secular theory which considers long-term aver- 
ages of forces between planets. This method was originally developed independently by 
Laplace and Lagrange, but see BG06a for a presentation that does not involve matrix 
manipulation. We compare the secular solutions to A^-body numerical models that solve 
the force equations directly, and are therefore arbitrarily accurate, for an sufficiently 
small timestep, using the MERCURY6 code (Chambers et aZ.1999). 

In general, secular theory predicts that e's and Atn's (the difference between two 
longitudes of periastron, w) oscillate and that a is constant. Therefore, conservation of 
angular momentum requires that as one planet's eccentricity drops, another must rise 
(as e increases, orbital angular momentum decreases). The type of oscillation depends on 
initial conditions. If Anj oscillates about 0, the system is experiencing "aligned libration." 
If Atu librates about tt, then the system is undergoing "anti-aligned libration." If Atn 
oscillates through 2tt then the apsides undergo "circulation". The boundaries between 
these qualitatively different types of behavior are known as "separatrices" . 

Recently it has been noted that many systems lie near an "apsidal separatrix" (Ford 
et aL2005; BG06a, BG06c). For systems of just two planets, the apsidal separatrix can 
only separate circulation and libration. This type of separatrix is a "libration-circulation 
separatrix" . An example of the libration-circulation separatrix is shown in the left panels 
of Fig. [H 

In systems of more than two planets, things get more complicated. In addition to 
the libration-circulation separatrix, the system may interact with different numbers of 
rotations of Atn through 360° during one eccentricity oscillation. The boundary be- 
tween interactions with different numbers of circulations in one eccentricity cycle is a 
"circulation-mode separatrix" , and an example is shown in the right panels of Fig. [TJ 
Note that the gray curve in the top right panel of this figure circulates once and librates 
once during one eccentricity cycle, whereas the black curve circulates twice. 

For an interaction to lie near a separatrix, the amplitude of eccentricity oscillations 
is generally two orders of magnitude or more. Since ^ e < 1 for bound planets, this 
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Figure 1. Examples of near-separatrix motion in planetary systems. Left Panels: A libration-cir- 
culation separatrix. Two possible evolutions of v And c and d (the middle and outer planet of the 
system) assuming the current orbits. The black points are the system from Butler et a/. (2006), 
the gray from Ford et o/.(2005). Although the best-fit orbits in these two cases are very similar, 
they result in qualitatively different types of evolution of Am (top): The older data predict 
aligned libration, whereas the updated data predict circulation. Note that the evolution of Cc is 
similar in both cases, and periodically reaches near-zero values (bottom). Right Panels: A cir- 
culation-mode separatrix. HD 69830 c and d evolve near the circulation-mode separatrix. The 
black data are from Lovis et aZ.(2006), and the gray data arc for a fictitious system in which the 
inner planet's, b's, eccentricity was changed from 0.1 to 0.15. In the first 2 x lO** years, the actual 
Azu undergoes 1 complete rotation through 360°, but in the fictitious system, Azu undergoes 
2 complete circulations (top). We again see that the middle planet's eccentricity periodically 
drops to near-zero values in both cases (bottom). 

means that at least one planet in near-separatrix interactions (both libration-circulation 
and circulation-mode) periodically is on a nearly circular orbit. The proximity to the 
separatrix can be parameterized as 

_ 2[min(-\/ X'^ + y^)] 

where x = 6162 sin(Acc7); y = 6162 cos(Aci7), and the subscripts min and max refer to 
the extreme values these variables attain over a complete secular cycle (BG06c). (Note 
that e is a fixed quantity even as e and Azu change periodically.) When e = the pair 
is on an apsidal separatrix, and one eccentricity periodically reaches zero. Note that in 
the case of mean motion resonant interactions, e may not necessarily reveal how close an 
interaction is to an apsidal separatrix (BG06c). 



(2.1) 



3. The Dynamical Stability Limit 

Observational uncertainties in the orbits of exoplanets often include regions of dynam- 
ical instability in which one or more planets would be ejected within 1 Myr, i.e. the 
system is dynamically "packed" (e.g. Barnes & Quinn 2001, 2004; Gozdziewski 2002; 
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Kiseleva-Eggleton 2002; Erdi et aL2004; BG06b, BG07b). This timescale appears to be 
long enough to identify nearly all unstable configurations (Barnes & Quinn 2004). Rep- 
resentative cases of these investigations are shown in Fig. [2l The shading in this figure 
correlates with the likelihood of stability: White regions are stable, black unstable, and 
grays contain a mix of stable and unstable orbits. In these cases we define stability in 
the "Lagrange" sense for ~ 10^ orbits: Escapes and exchanges of planets are forbidden. 

The preponderance of systems that lie near the stability limit naturally raises the ques- 
tion of where, exactly, in orbital element space, the stability boundary lies. Surprisingly, 
a general analytical expression for the location of the stability boundary is lacking, even 
though it has been the goal of considerable research (e.g. Szebehely & McKenzie 1981; 
Wisdom 1982; Marchal & Bozis 1982; Gladman 1993; Laskar 1990; Holman & Wiegert 
1999; Sosnitskii 1999; David et a/.2003; Cuntz et aL2007). These investigations all made 
significant progress toward defining a stability boundary, but none found it. 

Given that many planetary systems lie near the dynamical stability limit, we have 
considered how to quantify the proximity to stabihty (BG06b, BG07b). We will focus on 
one formulation by Marchal & Bozis (1982; see also Gladman 1993) that considers "Hill 
stability" . Hill stability only requires the ordering of the planets to remain constant for 
all time; the outer planet may escape to infinity. 

A coplanar, two-planet system, outside of a resonance is Hill stable if 

2^ ^ 1 , q4/3 "^i"^2 mim2(llmi -t-7TO2) , 



G^M^ ml^^{mi+m2)^^^ 3to3(™i + "^2)^ 

where M is the total mass of the system, mi is the mass of the more massive planet, 
7712 is the mass of the less massive planet, is the mass of the star, G is the gravita- 
tional constant, Af* — mim2 + rriims + 77127713, c is the total angular momentum of the 
system, and h is the energy (Marchal & Bozis 1982). Here c and h must be calculated in 
barycentric coordinates. If a given three-body system satisfies the inequality in Eq. (j3.ip . 
then the system is Hill stable. If this inequality is not satisfied, then the system may or 
may not be Hill stable. In this inequality, the left-hand side is a function of the orbits, 
but the right-hand side is only a function of the masses. This approach for identifying 
stable orbits is fundamentally different from other common techniques for determining 
stability which exploit resonance overlaps (Wisdom 1982; Quillen & Faber 2006), chaotic 
diffusion (Laskar 1990; Pepe et al.2007), fast Lyapunov indicators (Froeschle et a/.1997; 
Sandor et a/.2007), or periodic orbits (Voyatzis & Hadjidemetriou 2006; Hadjidemetriou 
2006). 

The ratio of the two sides of the inequality in Eq. p.ip (left-hand side over right-hand 
side) allows a quantification of a system's proximity to the Hill stability boundary. We 
define this ratio as /3. If /3 < 1 a system's Hill stability is unknown, if /3 > 1, Hill stability 
is guaranteed, if /3 = 1, the system is on the boundary. Note that for any system of two 
planets outside of resonance with fixed mass, the locus of Hill-stable orbits can be easily 
calculated. 

Fig. [2] includes contour lines of constant (3 for the four systems considered. For the 
bottom two cases (systems with no mean motion resonance), the (3 = 1 lines are near 
the transition from Lagrange stable orbits (white) to unstable (black). These two cases 
therefore show that when (3 >1, a system is not only Hill stable, it is also Lagrange stable. 

However, in the top two panels of Fig. [21 the (3 = 1 contour cuts through a swath of 
stability located at resonance. This qualitatively different behavior is no surprise because 
the criterion for Hill stability, Eq. (|3.ip . does not apply to systems in resonance. Note 
however that in the top left panel, outside of resonance, the /? contours do follow the 
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Figure 2. Lagrange stability boundary in relation to the Hill stability boundary for four exo- 
planetary systems. White regions represent bins in which all configurations were stable, black 
bins contained no stable configurations, darker shades of gray correspond to regions in which 
the fraction of stable simulations were smaller (e.g. Barnes & Quinn 2004). The curves represent 
contour lines of /3. Contour lines follow the shape of the Lagrange stability boundary, except in 
resonance, where the Lagrange stability region is extended to lower values of /3. 

Lagrange stability boundary, again at values near \. In the top right, we see that all (3 
values are less than 1, hence there are no non- resonant, Lagrange stable configurations 
in the region of orbital element space considered. 



4. Dynamical Properties of Exoplanets 

In this section we will tabulate the observed distributions of e and [3 and interpret 
these distributions in the context of models of planet formation. Observed extrasolar 
planets tend to lie near the apsidal separatrix. In the top panel of Fig. [3l we plot the 
observed distribution of e values (solid line) . These values were calculated through an A^- 
body integration that required energy conservation to be better than 1 part in 10** with 
a symplectic integration technique. Over 40% of known systems have e < 0.01, including 
two of the three giant-planet pairs in the Solar System. 
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We tabulate the dynamical properties of multi-planet systems in Table iQ, where we 
also list the sources of the initial conditions we used. In this table AM stands for "apsidal 
motion" and the possibilities are circulation (C), aligned libration (A), or anti-aligned 
libration (AA). The MMR column lists the resonance, if applicable. The proximities to 
the apsidal separatrix, e, are the values from the literature (BG06c), as are proximities to 
the Hill stability boundary, (3 (BG07b). The "Class" distinguishes orbits whose evolution 
is dominated by tidal (T), resonant (R) or secular (S) interactions. Table 1 includes the 
dynamical properties of the giant planets in our Solar System for comparison. 

Given the strong tendency for systems to be near the apsidal separatrix, it is natural to 
wonder how systems formed that way. The apsidal behavior is tied to the eccentricities. 
Considerable work by previous researchers has focused on reproducing the distribution 
of eccentricity values. The two most studied processes are "planet-planet scattering", 
i.e. gravitational encounters between planets (e.g. Rasio & Ford 1996; Weidenschilhng 
& Marzari 1996; Lin & Ida 1997; Ford et a/.2001; Marzari & Weidenschilhng 2002), and 
torquing by remnant protoplanctary disks that pumps up eccentricities (e.g. Artymowicz 
1992; Boss 2000; Papaloizou et a/.2001; Chiang & Murray 2002; Goldreich & Sari 2003; 
D'Angelo et a/.2006). Simulations of planet-planet scattering can produce eccentricity 
distributions that are similar to the observed distribution (Ford et aZ.2001), but models 
of disk torquing tend to produce eccentricities that are too low (D'Angelo et a/. 2006). 

Ford et aZ.(2005) were the first to consider how scattering may lead to near-separatrix 
behavior. They showed that the ejection of a hypothetical additional planet could produce 
the observed near-separatrix behavior among the planets orbiting v And. In their model 
three planets, initially on circular orbits, formed too close together to be stable. Scattering 
between two planets results in one of them being ejected. The third planet, orbiting at 
a safe distance from the scattering planets, remained on a circular orbit. The planet 
that remained bound to the star after the ejection received an impulsive kick and its 
eccentricity quickly jumped to a relatively large value. This sudden change created a 
new "initial condition" for the secular interaction. The planet that did not partake in 
the scattering then began a new secular evolution, but since apsidal behavior is periodic, 
that planet's eccentricity would return to zero and the motion was near-separatrix. 

However, Ford et aZ.only considered one such unstable case, and it was unclear how 
likely such a scenario was. The observed e statistics (solid line in the top panel of Fig. [3]) 
place a strong constraint on planet-planet scattering models. BG07a considered ^ 400 
hypothetical systems in order to provide a statistical test of the planet-planet scattering 
model's likelihood to reproduce the observed apsidal behavior of extrasolar planetary 
systems. In the main part of that experiment, we considered systems in which all orbital 
elements were initially the same except for the mean longitude of the middle planet (of 
a three-planet system) which was shifted by 1 degree for each successive case (i.e. 360 
simulations). Stability is most likely independent of this parameter (see Eg. 13. 1|) . there- 
fore there should be no correlation between initial mean longitude and e. The resulting 
distribution of e is shown in the top panel of Fig. [3] by the dashed line. This scattering 
(dashed line) produced one-tenth the fraction of near-separatrix (e < 0.01) cases as is 
observed (solid hue). 

The planet-planet scattering model had done a reasonable job of reproducing the ob- 
served eccentricity distribution (Ford et a/.2001), but it may be a challenge for it to 
reproduce the observed apsidal behavior (BG07a). BG07a noted that the path to small 
e values seemed quite constrained: The scattered planet needed to be ejected immedi- 



t see |http://www.lpl.arizona.edu/~rory/research/xsp/dynamics| for an up-to-date list of 
these properties 
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Figure 3. Top Panel: Distributions of e. Tlie solid line is tlie observed distribution, the dashed 
line is that predicted by a model involving the ejection of an additional Jupiter mass planet (Ford 
et aZ.2005; BG07a). Bottom Panel: Observed distribution of /3 in known two-planet systems. 
Note that we may only calculate f3 for two-planet systems. 



ately. This common scattering route to near-separatrix motion led BG07a to suggest 
that the probability of near-separatrix motion could be enhanced if the perturber is rel- 
atively small and/or on an initially eccentric orbit. Such a "rogue planet" may have a 
better chance of being removed from a planetary system after one encounter with a giant 
planet, but this model will require further testing. The actual distribution of e must be 
reproduced by any model of planet formation. 

Next we discuss the implications of the /3 distribution, that is, how close a system is 
to instability. We tabulate f3 (Table 1) values to produce another distribution function 
that describes the range of dynamical properties of planetary systems. The observed 
distribution of /3 is plotted in the bottom panel of Fig. [31 Most systems have (3^1, 
nearly all have /3 < 2, and the exception is HD 217107 {(3 > 7). This clustering of 
systems near the Hill stability limit suggests that many systems form near the edge of 
dynamical stability. 

In Fig. [3] and Table 1 we see that many systems in fact have /? < 1. These systems 
are all in mean motion resonance; the interaction is stabilized by the resonance. In fact. 
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HD 108874 is the only known resonance system to have (3 > 1. Therefore it seems that 
resonances tend to form such that (3 < 1. 



Table 1 Dynamical Properties of Multiple Planet Systems 



System 


Pair 


MMR 


AM 


e 





Class 


Reference 


47 UMa 


b-c 


- 


C 





1.025 


S 


Butler et a/. (2006) 


55 Cnc 


e-b 


- 


c 


0.067 


- 


T 


Butler et a/. (2006) 




b-c 


3:1 


c 


0.11 


- 


R 






c-d 


- 


c 


0.158 


- 


S 




GJ 317 


b-c 


4:1? 


? 


- 


0.98 


R? 


Johnson et al.{2007) 


GJ 876 


d-c 


- 







- 


T 


Butler et aL(2006) 




c-b 


2:1 


A 


0.34 


- 


R 




Gl 581 


b-c 


- 


C 


0.15 


- 


T 


Udry et al.{2007) 




c-d 


- 


C 


0.20 


- 


T 




HD 12661 


b-c 


- 


c 


0.003 


1.199 


S 


Butler et a/. (2006) 


HD 37124 


b-c 


- 


c 


0.009 


- 


S 


Butler et a/. (2006) 




c-d 


- 


A 


0.096 


- 


S 




HD 38529 


b-c 


- 


c 


0.44 


2.070 


s 


Butler et a/. (2006) 


HD 69830 


b-c 


- 


c 


0.095 


- 


T 


Lovis et a/. (2006) 




c-d 


- 


c 


0.04 


- 


S 




HD 74156'= 


b-c 


- 


c 


0.36 


1.542 


s 


Butler et a/. (2006) 


HD 73526 


b-c 


2:1 


AA 


0.006 


0.982 


R 


Butler et aZ.(2006) 


HD 82943 


b-c 


2:1 


C 


0.004 


0.946 


R 


Butler et aL(2006) 


HD 128311 


b-c 


2:1 


C 


0.091 


0.968 


R 


Butler et a/. (2006) 


HD 108874 


b-c 


4:1 


C/AA'^ 


0.2 


1.107 


R 


Butler et aZ.(2006) 


HD 155358 


b-c 


- 


AA 


0.21 


1.043 


S 


Cochran et al.{2007) 


HD 168443 


b-c 


- 


C 


0.22 


1.939 


S 


Butler et a/. (2006) 


HD 169830 


b-c 


- 


C 


0.33 


1.280 


s 


Butler et d.(2006) 


HD 190360 


c-b 


- 


c 


0.38 


1.701 


T 


Butler et aL(2006) 


HD 202206 


b-c 


5:1 


c 


0.096 


0.883 


R 


Butler et a/. (2006) 


HD 217107 


b-c 


- 


c 


0.46 


7.191 


T 


Butler et aL(2006) 


Hip 14810 


b-c 


- 


AA 


0.05 


1.202 


T 


Wright et aL(2007) 


SS 


J-S 




c 


0.19 




S 


JPL 




S-U 




c 


0.006 




S 






U-N 




c 


0.004 




s 




/U Ara 


c-d 




c 


0.002 




T 


Pepe et aL(2007) 




d-b 


2:1 


c 


0.003 




R 






b-e 




c 


0.13 




s 




V And 


b-c 




c 


1.8 X 10"" 




T 


Butler et aL(2006) 




c-d 




c 


2.8 X 10"^ 




S 





" The current eccentricity of one planet is 0, placing the pair on an apsidal separatrix. 

^ The values of mean longitude and longitude of pcriastron arc unknown. 
These values do not incorporate the new planet Bean et aL(2007) 
This pair alternates between circulation and anti-aligned libration. 

One noteworthy system in this context is GJ 317 (Johnson et aL2007), with a (3 value 
of 0.98 (see Table 1). This value is only permissible if the system is in a mean motion 
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resonance, and, indeed, the ratio of the periods is near a whole number (4.02). UsuaUy a 
mean motion resonance is identified by calculating the "resonant argument" , which de- 
pends on mean longitude and longitude of periastron (see, e.g. Murray & Dermott 1999). 
However, in the case of GJ 317, these two angles are unknown, and, hence, the resonant 
argument is also unknown. Nonetheless, the system must be in resonance (barring sub- 
stantial revisions in the masses, semi-major axes and eccentricities of the planets) based 
on our stability analysis. Therefore consideration of proximity to instability in terms of 
the Hill stability boundary may provide an alternative method for identifying mean mo- 
tion resonances. As the orbital parameters of the two planets in GJ 317 are refined, we 
predict that they will, in fact, be found to be in resonance. 

5. Conclusions 

We have described here new approaches for parameterizing the dynamical interactions 
of extrasolar planet systems. Many systems interact such that they arc near the stability 
limit and the apsidal scparatrix. Most mean-motion resonance interactions would be 
unstable if not for the resonance, a feature which may be used to identify resonant 
interactions (e.g. GJ 317 b and c). The distribution of e shows that the scattering of 
approximately equal mass planets is unlikely to produce the observed distribution of 
apsidal behavior. The distribution of /? values shows that many systems formed near the 
limit of dynamical stability. 

These results demonstrate the benefits of the consideration of the dynamical interac- 
tions of multiple planet systems, which may constrain models of planet formation. For 
example, as we discussed above, the proximity of so many systems to an apsidal sepa- 
ratrix (e = 0) suggests that "rogue planets" may have played a major role in scattering 
planets to higher eccentricities. The proximity of so many systems to the stability limit 
(/? = 1) suggests many systems form in densely packed configurations. Moreover, con- 
sideration of planetary systems in these terms has revealed that our Solar System shares 
dynamical traits with the known multiple planet systc;ins. Perhaps constraining planet 
formation models through dynamical properties, rather than observed orbital elements, 
will lead to a universal model of planet formation. 

About half of planetary systems are multiple (Wright et aZ.2007), predictions of addi- 
tional companions are being borne out (Barnes & Raymond 2004; Raymond & Barnes 
2005; Bean et a/.2008), and the current distribution of planet masses suggest there will 
be many planets with a mass equal to that of Saturn or less (Marcy et a/.2005), i.e. below 
current detection limits. These three observations imply many multiple planet systems 
will be detected in the future. Hence characterizing extrasolar planet interactions will be 
a critical aspect of the study of planet formation for the foreseeable future. 
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